rm(list = ls())
getwd()
setwd("C:/Users/griff/OneDrive/Desktop")

require(pscl)
require(MASS)
require(boot)


data <- read.csv(file="Liberal Constraints/JTPV/Data/SQ_attacks_during_insurgency.csv", na.string = "")
str(data)
data$IGO_Context <- as.numeric(as.character(data$IGO_Context))
data$loggdp <- as.numeric(as.character(data$loggdp))
data$polity <- as.numeric(as.character(data$polity))


n1 <- zeroinfl(SQ_attacks ~ Democracy + IGO_Context + softpta +
                 territorial_control + loggdp + past_incidents,
               data=data, dist = "negbin")
summary(n1)

n2 <- zeroinfl(SQ_attacks ~ Democracy + IGO_Context +
                 territorial_control + loggdp + past_incidents,
               data=data, dist="negbin")
summary(n2)

n3 <- zeroinfl(SQ_attacks ~ Democracy + softpta +
                 data=data, dist="negbin")
summary(n3)

n4 <- zeroinfl(SQ_attacks ~ polity + IGO_Context + softpta +
                 territorial_control + loggdp + past_incidents,
               data=data, dist="negbin")
summary(n4)

n5 <- zeroinfl(SQ_attacks ~ e_boix_regime + 
                 territorial_control + loggdp + past_incidents,
               data=data, dist="negbin")

summary(n5)


#now comparing the current model with the null model without predictors
#using a chi-squared test on the difference of log likelihoods

m0 <- update(n1, . ~ 1)
pchisq(2 * (logLik(n1) - logLik(m0)), df = 3, lower.tail=FALSE)
#Log Lik 8.873981e-26 (df=12) we can tell that our overall model is statistically significant

#NOW have to tell if the zero-inflated model is an improvement over a regular negative binomial

m2 <- glm.nb(SQ_attacks ~ e_boix_regime + IGOc + softpta +
               territorial_control + loggdp + past_incidents, data = insurgencyset)
summary(m2)

vuong(n1, m2)
#RAW model1 > model2 pvalue: 0.031532
#AIC-corrected model1 > model2 pvalue: 0.067113
#BIC-corrected model1 > model2 pvalue 0.264539
